##I Rimbalzi Nel Basket
library(here)
library(corrplot)
library(nortest)
library(lmtest)
library(car)
library(plotly)
library(heatmaply)
library(ggheatmap)
library(ggplot2)
library(viridis)
library(glmnet)
library(highcharter)
library(reshape2)
library(RColorBrewer)
library(tidyr)
library(dplyr)
knitr::opts_knit$set(root.dir = here("0_Materiale"))
NBA moderna (1976-2011): VARIABILE DIPENDENTE: numero di vittorie in stagione COVARIATE: tutte le altre (o uno specifico insieme di queste, in base all’obiettivo di analisi) Considerare solo le squadre che hanno giocato 82 partite (dataset$games==82)
dataset <- read.delim("basketball_teams.txt") # andiamo a leggere il database fornito
FIRST <- 1976 # primo anno del range da considerare per lo studio
LAST <- 2011 # ultimo anno del range da considerare per lo studio
df <- dataset [dataset$lgID=="NBA" & dataset$year >= FIRST & dataset$year <= LAST & dataset$games==82,]
dataset$lgID <- as.factor(dataset$lgID) # perchè mi permettono di poter generare variabili dummy
summary(df)
## year lgID tmID franchID
## Min. :1976 Length:841 Length:841 Length:841
## 1st Qu.:1985 Class :character Class :character Class :character
## Median :1993 Mode :character Mode :character Mode :character
## Mean :1993
## 3rd Qu.:2001
## Max. :2008
## confID divID rank confRank
## Length:841 Length:841 Min. :0.000 Min. : 1.000
## Class :character Class :character 1st Qu.:2.000 1st Qu.: 4.000
## Mode :character Mode :character Median :3.000 Median : 7.000
## Mean :3.565 Mean : 7.164
## 3rd Qu.:5.000 3rd Qu.:10.000
## Max. :8.000 Max. :15.000
## playoff name o_fgm o_fga
## Length:841 Length:841 Min. :2565 Min. :5972
## Class :character Class :character 1st Qu.:2981 1st Qu.:6592
## Mode :character Mode :character Median :3220 Median :6903
## Mean :3239 Mean :6941
## 3rd Qu.:3489 3rd Qu.:7253
## Max. :3980 Max. :8868
## o_ftm o_fta o_3pm o_3pa o_oreb
## Min. :1189 Min. :1475 Min. : 0.0 Min. : 0.0 Min. : 720
## 1st Qu.:1523 1st Qu.:2039 1st Qu.: 78.0 1st Qu.: 293.0 1st Qu.: 974
## Median :1649 Median :2201 Median :283.0 Median : 814.0 Median :1083
## Mean :1666 Mean :2212 Mean :279.1 Mean : 804.6 Mean :1086
## 3rd Qu.:1797 3rd Qu.:2371 3rd Qu.:445.0 3rd Qu.:1267.0 3rd Qu.:1191
## Max. :2388 Max. :3051 Max. :837.0 Max. :2283.0 Max. :1520
## o_dreb o_reb o_asts o_pf o_stl
## Min. :2044 Min. :2922 Min. :1422 Min. :1476 Min. : 455.0
## 1st Qu.:2348 1st Qu.:3381 1st Qu.:1760 1st Qu.:1777 1st Qu.: 613.0
## Median :2433 Median :3506 Median :1934 Median :1900 Median : 671.0
## Mean :2434 Mean :3520 Mean :1934 Mean :1909 Mean : 681.3
## 3rd Qu.:2527 3rd Qu.:3647 3rd Qu.:2094 3rd Qu.:2033 3rd Qu.: 746.0
## Max. :2966 Max. :4216 Max. :2575 Max. :2470 Max. :1059.0
## o_to o_blk o_pts d_fgm d_fga
## Min. : 910 Min. :204.0 Min. : 6901 Min. :2488 Min. :5638
## 1st Qu.:1207 1st Qu.:360.0 1st Qu.: 7958 1st Qu.:2978 1st Qu.:6593
## Median :1311 Median :411.0 Median : 8404 Median :3243 Median :6911
## Mean :1338 Mean :421.8 Mean : 8423 Mean :3239 Mean :6941
## 3rd Qu.:1443 3rd Qu.:473.0 3rd Qu.: 8879 3rd Qu.:3493 3rd Qu.:7268
## Max. :2011 Max. :716.0 Max. :10371 Max. :4265 Max. :8142
## d_ftm d_fta d_3pm d_3pa d_oreb
## Min. :1217 Min. : 0.0 Min. : 0.0 Min. :1579 Min. : 745
## 1st Qu.:1514 1st Qu.: 82.0 1st Qu.: 282.0 1st Qu.:2033 1st Qu.: 986
## Median :1648 Median :280.0 Median : 836.0 Median :2203 Median :1095
## Mean :1666 Mean :279.1 Mean : 804.6 Mean :2212 Mean :1086
## 3rd Qu.:1808 3rd Qu.:446.0 3rd Qu.:1251.0 3rd Qu.:2395 3rd Qu.:1177
## Max. :2377 Max. :683.0 Max. :1768.0 Max. :3071 Max. :1495
## d_dreb d_reb d_asts d_pf d_stl
## Min. :2012 Min. :2976 Min. :1336 Min. :1434 Min. :461.0
## 1st Qu.:2326 1st Qu.:3378 1st Qu.:1778 1st Qu.:1788 1st Qu.:623.0
## Median :2431 Median :3497 Median :1939 Median :1900 Median :677.0
## Mean :2434 Mean :3520 Mean :1934 Mean :1909 Mean :681.3
## 3rd Qu.:2529 3rd Qu.:3653 3rd Qu.:2092 3rd Qu.:2020 3rd Qu.:734.0
## Max. :3067 Max. :4309 Max. :2537 Max. :2453 Max. :955.0
## d_to d_blk d_pts o_tmRebound d_tmRebound
## Min. : 949 Min. :264.0 Min. : 6909 Min. :0 Min. :0
## 1st Qu.:1208 1st Qu.:380.0 1st Qu.: 7968 1st Qu.:0 1st Qu.:0
## Median :1304 Median :419.0 Median : 8453 Median :0 Median :0
## Mean :1338 Mean :421.8 Mean : 8423 Mean :0 Mean :0
## 3rd Qu.:1444 3rd Qu.:460.0 3rd Qu.: 8841 3rd Qu.:0 3rd Qu.:0
## Max. :1980 Max. :654.0 Max. :10723 Max. :0 Max. :0
## homeWon homeLost awayWon awayLost neutWon
## Min. : 6.00 Min. : 1.00 Min. : 1.00 Min. : 8.00 Min. :0
## 1st Qu.:21.00 1st Qu.:10.00 1st Qu.:10.00 1st Qu.:21.00 1st Qu.:0
## Median :26.00 Median :15.00 Median :15.00 Median :26.00 Median :0
## Mean :25.63 Mean :15.37 Mean :15.37 Mean :25.63 Mean :0
## 3rd Qu.:31.00 3rd Qu.:20.00 3rd Qu.:20.00 3rd Qu.:31.00 3rd Qu.:0
## Max. :40.00 Max. :35.00 Max. :33.00 Max. :40.00 Max. :0
## neutLoss confWon confLoss divWon divLoss
## Min. :0 Min. : 5.00 Min. : 7.00 Min. : 1.00 Min. : 1.00
## 1st Qu.:0 1st Qu.:20.00 1st Qu.:20.00 1st Qu.: 8.00 1st Qu.: 9.00
## Median :0 Median :27.00 Median :26.00 Median :12.00 Median :12.00
## Mean :0 Mean :26.91 Mean :26.91 Mean :12.27 Mean :12.27
## 3rd Qu.:0 3rd Qu.:34.00 3rd Qu.:33.00 3rd Qu.:16.00 3rd Qu.:16.00
## Max. :0 Max. :48.00 Max. :52.00 Max. :25.00 Max. :27.00
## pace won lost games min
## Min. : 0.00 Min. :11 Min. :10 Min. :82 Min. :19680
## 1st Qu.: 0.00 1st Qu.:31 1st Qu.:32 1st Qu.:82 1st Qu.:19780
## Median : 0.00 Median :42 Median :40 Median :82 Median :19805
## Mean : 6.71 Mean :41 Mean :41 Mean :82 Mean :19817
## 3rd Qu.: 0.00 3rd Qu.:50 3rd Qu.:51 3rd Qu.:82 3rd Qu.:19855
## Max. :102.00 Max. :72 Max. :71 Max. :82 Max. :20080
## arena attendance bbtmID
## Length:841 Min. : 0 Length:841
## Class :character 1st Qu.:32767 Class :character
## Mode :character Median :32767 Mode :character
## Mean :32728
## 3rd Qu.:32767
## Max. :32767
#STUDIO DATI RELATIVI AI RIMBALZI
# o_oreb = Rimbalzi ottenuti in attacco
# o_dreb = Rimbalzi subiti in attacco
# o_reb = totale rimbalzi in attacco
# d_oreb = Rimbalzi subiti in difesa
# d_dreb = Rimbalzi ottenuti in difesa
# d_reb = totale rimbalzi in difesa
df_reb <- subset(df, select = c("o_oreb", "o_dreb", "o_reb", "d_oreb", "d_dreb", "d_reb", "won"))
par(mfrow = c(2, 3)) # Imposta il layout a 2 righe e 3 colonne
for (variables in 1:(dim(df_reb)[2]-1)){
thisvar = df_reb[,variables]
d <- density(thisvar)
xmin <- floor(min(thisvar))
xmax <- ceiling(max(thisvar))
# Crea il plot della densità con stile accattivante
plot(d, main = names(df_reb)[variables], xlab = "", col = "skyblue", lwd = 1, xlim = c(xmin, xmax), ylim = c(0, max(d$y)*1.1))
# Aggiungi la distribuzione normale teorica ideale in rosso
x <- seq(xmin, xmax, length = 100)
lines(x, dnorm(x, mean = mean(thisvar), sd = sd(thisvar)), col = "red", lwd = 1)
# Aggiungi griglia per migliorare la leggibilità
grid()
}
# Aggiungi titolo
title("Density plots with Normal Distribution", line = -17, cex.main = 2, outer = TRUE)
#GRAFICO INTUITIVO DA TESTARE
# Imposta il layout a 2 righe e 3 colonne
#par(mfrow = c(2, 3))
# Lista per salvare i plot interattivi
#plot_list <- list()
#for (variables in 1:(dim(df_reb)[2])) {
# thisvar <- df_reb[, variables]
# d <- density(thisvar)
# xmin <- floor(min(thisvar))
# xmax <- ceiling(max(thisvar))
# Crea il plot della densità con plot_ly
# p <- plot_ly(x = d$x, y = d$y, type = "scatter", mode = "lines", line = list(color = "skyblue", width = 1),
# name = names(df_reb)[variables]) %>%
# add_trace(x = x, y = dnorm(x, mean = mean(thisvar), sd = sd(thisvar)), type = "scatter", mode = "lines",
# line = list(color = "red", width = 1), name = "Normal Distribution") %>%
# layout(title = names(df_reb)[variables])
# Aggiungi il plot alla lista
# plot_list[[variables]] <- p
#}
# Aggiungi titolo
#title("Density plots with Normal Distribution", line = -17, cex.main = 2, outer = TRUE)
# Visualizza i plot interattivi
#subplot(plot_list)
#verifica della presenza zeri nel dataset
print(paste("Percentage non-zero o_oreb: ",round(length(which(df_reb$o_oreb>0)) /dim(df_reb)[1]*100,2)))
## [1] "Percentage non-zero o_oreb: 100"
print(paste("Percentage non-zero o_dreb: ",round(length(which(df_reb$o_dreb>0)) /dim(df_reb)[1]*100,2)))
## [1] "Percentage non-zero o_dreb: 100"
print(paste("Percentage non-zero o_reb: ",round(length(which(df_reb$o_reb>0)) /dim(df_reb)[1]*100,2)))
## [1] "Percentage non-zero o_reb: 100"
print(paste("Percentage non-zero d_oreb: ",round(length(which(df_reb$d_oreb>0)) /dim(df_reb)[1]*100,2)))
## [1] "Percentage non-zero d_oreb: 100"
print(paste("Percentage non-zero d_dreb: ",round(length(which(df_reb$d_dreb>0)) /dim(df_reb)[1]*100,2)))
## [1] "Percentage non-zero d_dreb: 100"
print(paste("Percentage non-zero d_reb: ",round(length(which(df_reb$d_reb>0)) /dim(df_reb)[1]*100,2)))
## [1] "Percentage non-zero d_reb: 100"
#STUDIO SULLE VITTORIE
M <- cor(as.matrix(df[, c(11:25, 54)])) # correlation matrix
corrplot(M, method="color", outline = TRUE,type="lower",order = "hclust",
tl.col="black", tl.srt=45, diag=FALSE,tl.cex = 1,mar=c(0,0,3,0),
title="Correlation Matrix between Predictor and Outcome variables")
#ISTOGRAMMA
# Crea un istogramma di base
hist(df$won, col = "skyblue", border = "white", main = "Distribuzione delle Vittorie", xlab = "Numero di Vittorie", ylab = "Frequenza")
# Aggiungi una griglia di sfondo
grid()
# Aggiungi una linea di riferimento
abline(v = mean(df$won), col = "red", lwd = 2)
# Aggiungi una legenda
legend("topright", legend = c("Media"), col = c("red"), lwd = 2)
#ISTOGRAMMA INTERATTIVO DA TESTARE 1° versione
# Crea un istogramma con plot_ly
histogram <- plot_ly(df, x = ~won, type = "histogram", marker = list(color = "skyblue", line = list(color = "white", width = 0.5))) %>%
layout(title = "Distribuzione delle Vittorie",
xaxis = list(title = "Numero di Vittorie"),
yaxis = list(title = "Frequenza"),
barmode = "overlay") %>%
add_trace(x = ~mean(won), type = "scatter", mode = "lines", line = list(color = "red", width = 2), name = "Media")
# Visualizza l'istogramma interattivo
histogram
## A marker object has been specified, but markers is not in the mode
## Adding markers to the mode...
#ISTOGRAMMA INTERATTIVO DA TESTARE 2° versione
# Crea un istogramma con plot_ly
histogram <- plot_ly(df, x = ~won, type = "histogram",
marker = list(color = "skyblue", line = list(color = "white", width = 0.5)),
opacity = 0.7) %>%
layout(title = list(text = "Distribuzione delle Vittorie", x = 0.5),
xaxis = list(title = "Numero di Vittorie", zeroline = FALSE),
yaxis = list(title = "Frequenza", zeroline = FALSE),
barmode = "overlay") %>%
add_trace(x = ~mean(won), type = "scatter", mode = "lines",
line = list(color = "red", width = 2), name = "Media") %>%
add_trace(x = ~mean(won), type = "scatter", mode = "markers",
marker = list(color = "red", size = 8), showlegend = FALSE) %>%
add_annotations(text = sprintf("Media: %.2f", mean(df$won)), x = mean(df$won), y = 0,
arrowhead = 2, arrowcolor = "red", arrowsize = 1.5, arrowwidth = 2,
ax = 0, ay = -40, xshift = 10, yshift = -10)
# Visualizza l'istogramma interattivo
histogram
## A marker object has been specified, but markers is not in the mode
## Adding markers to the mode...
# PLOT DENSITA
density_plot <- density(df$won)
# Plot di base
plot(density_plot, main = "Distribuzione di Densità delle Vittorie", col = "skyblue", lwd = 2, ylim = c(0, 0.07), xlim = c(0, max(df$won)))
# Aggiungi titoli ed etichette degli assi
title(main = "Distribuzione di Densità delle Vittorie", xlab = "Numero di Vittorie", ylab = "Densità")
# Aggiungi una griglia di sfondo
grid()
# Aggiungi una linea di riferimento per la media
abline(v = mean(df$won), col = "red", lwd = 2)
# Aggiungi una legenda
legend("topright", legend = c("Media"), col = c("red"), lwd = 2)
# PLOT DENSITA INTERATTIVO DA TESTARE
# Calcola la densità delle vittorie
density_plot <- density(df$won)
# Crea il plot di densità con plot_ly
density_interactive <- plot_ly(x = density_plot$x, y = density_plot$y, type = "scatter", mode = "lines",
line = list(color = "skyblue", width = 2)) %>%
layout(title = list(text = "Distribuzione di Densità delle Vittorie", x = 0.5),
xaxis = list(title = "Numero di Vittorie"),
yaxis = list(title = "Densità"),
showlegend = FALSE) %>%
add_trace(x = c(mean(df$won), mean(df$won)), y = c(0, max(density_plot$y)),
type = "scatter", mode = "lines", line = list(color = "red", width = 2),
name = "Media") %>%
add_trace(x = mean(df$won), y = max(density_plot$y), type = "scatter", mode = "markers",
marker = list(color = "red", size = 8)) %>%
add_annotations(text = sprintf("Media: %.2f", mean(df$won)), x = mean(df$won), y = max(density_plot$y) * 1.05,
arrowhead = 2, arrowcolor = "red", arrowsize = 1.5, arrowwidth = 2,
ax = 0, ay = -40, xshift = 10, yshift = -10)
# Visualizza il plot di densità interattivo
density_interactive
## A line object has been specified, but lines is not in the mode
## Adding lines to the mode...
#CORRPLOT
# # Calcola la matrice di correlazione
# M <- cor(as.matrix(df[, c(11:25, 54)]))
#
# # Crea il grafico della matrice di correlazione
# corrplot(M, method = "color", outline = TRUE, type = "lower", order = "hclust",
# tl.col = "black", tl.srt = 45, diag = FALSE, tl.cex = 0.7, mar = c(0, 0, 3, 0),
# title = "Correlation Matrix between Predictor and Outcome variables",
# addCoef.col = "black", number.cex = 0.5, tl.offset = 0.5)
#
# # Aggiungi una legenda
# legend("bottomleft", legend = c("Positive Correlation", "Negative Correlation"), fill = c("white", "black"), border = NA, bty = "n", title = "Correlation")
#
# CORRPLOT
data_subset <- df[, c(11:40, 54)]
cor_matrix <- cor(data_subset, use = "complete.obs")
cor_data <- melt(cor_matrix)
names(cor_data) <- c("Variable1", "Variable2", "Correlation")
# Create the heatmap
p <- plot_ly(data = cor_data, x = ~Variable1, y = ~Variable2, z = ~Correlation,
type = "heatmap",
colors = viridis(n = 1024, option = "magma"),
hoverinfo = "x+y+z") %>%
layout(title = 'Correlation Matrix Heatmap',
xaxis = list(tickangle = 45),
yaxis = list(tickangle = 45),
autosize = TRUE)
# cor_values <- round(as.matrix(cor_matrix), 2) # Round for readability
# for (i in seq_len(nrow(cor_matrix))) {
# for (j in seq_len(ncol(cor_matrix))) {
# p <- p %>% add_annotations(
# x = rownames(cor_matrix)[i],
# y = colnames(cor_matrix)[j],
# text = as.character(cor_values[i, j]),
# showarrow = FALSE,
# font = list(color = ifelse(cor_values[i, j] < 0.5, "white", "black"))
# )
# }
# }
(p)
#BOXPLOT
aggregated_data <- df %>%
select(tmID, year, won) %>%
group_by(tmID, year) %>%
summarize(
won = mean(won, na.rm = TRUE),
)
## `summarise()` has grouped output by 'tmID'. You can override using the
## `.groups` argument.
reshaped_data <- aggregated_data %>%
pivot_longer(cols = c(won), names_to = "stat_type", values_to = "stat") %>%
unite("new_col", tmID, stat_type, sep = "_") %>%
pivot_wider(names_from = new_col, values_from = "stat")
fig <- plot_ly()
for (i in 2:ncol(reshaped_data)) {
current_column_data <- reshaped_data[[i]]
fig <- fig %>% add_trace(y = current_column_data, name = colnames(reshaped_data)[i], type = "box")
}
(fig)
## Warning: Ignoring 30 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 19 observations
## Warning: Ignoring 4 observations
## Warning: Ignoring 23 observations
## Warning: Ignoring 8 observations
## Warning: Ignoring 24 observations
## Warning: Ignoring 12 observations
## Warning: Ignoring 13 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 29 observations
## Warning: Ignoring 30 observations
## Warning: Ignoring 31 observations
## Warning: Ignoring 31 observations
## Warning: Ignoring 13 observations
## Warning: Ignoring 9 observations
## Warning: Ignoring 26 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 19 observations
## Warning: Ignoring 3 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 21 observations
## Warning: Ignoring 11 observations
df$reb <- df$o_reb + df$d_reb
ad.test(df$reb)
##
## Anderson-Darling normality test
##
## data: df$reb
## A = 3.6997, p-value = 3.1e-09
Con un livello di significatività (\(\alpha\)) di 0.01 e un p-value molto piccolo (3.1e-09) ottenuto dal test di normalità di Anderson-Darling per i dati della variabile df$reb, puoi concludere che hai sufficiente evidenza statistica per respingere lipotesi nulla che i dati seguono una distribuzione normale.Con il tuo livello di significatività del 0.01 e il p-value molto piccolo (3.1e-09), il p-value è inferiore al livello di significatività, quindi respingeresti lipotesi nulla. Questo suggerisce che i dati nella variabile df$reb non seguono una distribuzione normale al livello di significatività del 0.01. In termini più pratici, hai abbastanza evidenza statistica per concludere che la variabile df$reb non segue una distribuzione normale basandoti sui risultati del test di Anderson-Darling.
ks.test(df$reb, "pnorm")
## Warning in ks.test.default(df$reb, "pnorm"): i legami non dovrebbero essere
## presenti per il test di Kolmogorov-Smirnov
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: df$reb
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided
Il risultato che hai ottenuto riguarda il test di Kolmogorov-Smirnov a campione singolo sui dati contenuti nella variabile df$reb. Il test KS confronta la distribuzione empirica dei tuoi dati con una distribuzione teorica (spesso una distribuzione uniforme). In breve, il risultato suggerisce che i tuoi dati non seguono la distribuzione teorica presunta, e cè un elevata probabilità che la differenza osservata sia statisticamente significativa.
sf.test(df$reb)
##
## Shapiro-Francia normality test
##
## data: df$reb
## W = 0.98016, p-value = 1.758e-08
In sintesi, il risultato del test di Shapiro-Francia indica che i tuoi dati nella variabile df$reb non seguono una distribuzione normale. Questo è supportato dal valore basso del p-value, il quale suggerisce che la differenza tra la distribuzione dei tuoi dati e una distribuzione normale è statisticamente significativa.
\(\text{Formula1} = \frac{\text{Rimbalzi offensivi in attacco}}{\text{Tiri sbagliati su azione}}\) Rappresenta la capacità della squadra di ripossesso della palla dopo un tiro che non va a canestro e colpisce il tabellone.
\(\text{Formula2} = \frac{\text{Rimbalzi difensivi in difesa presi}}{\text{Tiri sbagliati su azione degli avversari}}\) Rappresenta la capacità della squadra di impossessarsi della palla dopo un tiro sbagliato della squadra avversaria che colpisce il tabellone, che troviamo un buon stimatore della capacità di contropiede della squadra.
\(\text{Formula3} = \frac{\text{Palle riprese in attacco} + 1.5 \times \text{Palle riprese in difesa}}{\text{Palle perse in attacco} + 2 \times \text{Rimbalzi subiti in difesa}}\) Rappresenta il rapporto tra le palle riprese nei rimbalzi (sia offensivi che difensivi) rispetto alle palle perse nei rimbalzi (sia offensivi che difensivi). I coefficienti sono stati scelti in base a ciò che riteniamo più importante in una partita, ossia la difesa del proprio canestro.
\(\text{Formula4} = (\text{Palle riprese in attacco - Palle perse in attacco}) + 1.5*(\text{Palle riprese in difesa - Palle perse in difesa})\) Cresce all’aumentare dei rimbalzi ottenuti e diminuisce all’aumentare dei rimbalzi subiti, considerando anche un coefficiente che da particolare importanza alla difesa.
\(\text{Formula5} = \frac{(\frac{\text{Rimbalzi subiti in difesa}}{\text{Palle perse in difesa}})}{(\frac{\text{Rimbalzi subiti in attacco}}{\text{Palle perse in attacco}})}\) Mostra quanto siano influenti i rimbalzi nel rapporto tra le palle perse dalla squadra e le palle perse dagli avversari.
# o_oreb = Rimbalzi ottenuti in attacco
# o_dreb = Rimbalzi subiti in attacco
# o_reb = totale rimbalzi in attacco
# d_oreb = Rimbalzi subiti in difesa
# d_dreb = Rimbalzi ottenuti in difesa
# d_reb = totale rimbalzi in difesa
df$f1 <- (df$o_oreb)/(df$o_fga-df$o_fgm)
df$f2 <- (df$d_dreb)/(df$d_fga-df$d_fgm)
df$f3 <- (df$o_oreb + 1.5 * df$d_dreb)/(df$o_dreb + 2 * df$d_oreb)
df$f4 <- (df$o_oreb - df$o_dreb) + 1.5 * (df$d_dreb - df$d_oreb)
df$f5 <- (df$d_oreb / df$d_to) / (df$o_dreb / df$o_to)
linMod <- lm(won ~ f1 + f2 + f3 + f4 + f5, data = df)
summary (linMod)
##
## Call:
## lm(formula = won ~ f1 + f2 + f3 + f4 + f5, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.624 -3.796 0.044 3.243 17.608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 395.45317 18.36490 21.533 < 2e-16 ***
## f1 83.35460 10.45130 7.976 4.99e-15 ***
## f2 -59.76252 13.37858 -4.467 9.02e-06 ***
## f3 -273.63124 15.89419 -17.216 < 2e-16 ***
## f4 0.03391 0.00415 8.171 1.13e-15 ***
## f5 -178.48989 3.87823 -46.024 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.985 on 835 degrees of freedom
## Multiple R-squared: 0.8401, Adjusted R-squared: 0.8392
## F-statistic: 877.6 on 5 and 835 DF, p-value: < 2.2e-16
# Installa e carica le librerie necessarie
if (!requireNamespace("ggplot2", quietly = TRUE)) {
install.packages("ggplot2")
}
if (!requireNamespace("gridExtra", quietly = TRUE)) {
install.packages("gridExtra")
}
library(ggplot2)
library(gridExtra)
## Warning: il pacchetto 'gridExtra' è stato creato con R versione 4.3.2
##
## Caricamento pacchetto: 'gridExtra'
## Il seguente oggetto è mascherato da 'package:dplyr':
##
## combine
# Funzione per creare il plot di Residuals vs Fitted (Residui vs. Valori Previsti)
res_vs_fitted <- ggplot(linMod, aes(.fitted, .resid)) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Residuals vs Fitted", x = "Fitted Values", y = "Residuals") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Funzione per creare il plot QQ dei Residui
qq_plot <- ggplot(linMod, aes(sample = .stdresid)) +
stat_qq(color = "skyblue", size = 2) +
stat_qq_line(color = "red", size = 1) +
labs(title = "QQ Plot dei Residui", x = "Quantiles teorici", y = "Quantiles osservati") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot Scale-Location
scale_location_plot <- ggplot(linMod, aes(.fitted, sqrt(abs(.stdresid)))) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Scale-Location Plot", x = "Fitted Values", y = "sqrt(|Standardized Residuals|)") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot Residuals vs Leverage (Residui vs. Leverage)
res_vs_leverage <- ggplot(linMod, aes(.hat, .stdresid)) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Residuals vs Leverage", x = "Leverage", y = "Standardized Residuals") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Visualizza i quattro plot
grid.arrange(res_vs_fitted, qq_plot, scale_location_plot, res_vs_leverage, ncol = 2)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#PLOT INTERATTIVI DA TESTARE
# Installa e carica le librerie necessarie
if (!requireNamespace("ggplot2", quietly = TRUE)) {
install.packages("ggplot2")
}
if (!requireNamespace("gridExtra", quietly = TRUE)) {
install.packages("gridExtra")
}
if (!requireNamespace("plotly", quietly = TRUE)) {
install.packages("plotly")
}
library(ggplot2)
library(gridExtra)
library(plotly)
# Funzione per creare il plot di Residuals vs Fitted (Residui vs. Valori Previsti)
res_vs_fitted <- ggplot(linMod, aes(.fitted, .resid)) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Residuals vs Fitted", x = "Fitted Values", y = "Residuals") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot QQ dei Residui
qq_plot <- ggplot(linMod, aes(sample = .stdresid)) +
stat_qq(color = "skyblue", size = 2) +
stat_qq_line(color = "red", size = 1) +
labs(title = "QQ Plot dei Residui", x = "Quantiles teorici", y = "Quantiles osservati") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot Scale-Location
scale_location_plot <- ggplot(linMod, aes(.fitted, sqrt(abs(.stdresid)))) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Scale-Location Plot", x = "Fitted Values", y = "sqrt(|Standardized Residuals|)") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot Residuals vs Leverage (Residui vs. Leverage)
res_vs_leverage <- ggplot(linMod, aes(.hat, .stdresid)) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Residuals vs Leverage", x = "Leverage", y = "Standardized Residuals") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Converti i singoli grafici in oggetti plotly
plotly_res_vs_fitted <- ggplotly(res_vs_fitted)
## `geom_smooth()` using formula = 'y ~ x'
plotly_qq_plot <- ggplotly(qq_plot)
plotly_scale_location_plot <- ggplotly(scale_location_plot)
## `geom_smooth()` using formula = 'y ~ x'
plotly_res_vs_leverage <- ggplotly(res_vs_leverage)
## `geom_smooth()` using formula = 'y ~ x'
# Visualizza i quattro grafici interattivi
subplot(plotly_res_vs_fitted, plotly_qq_plot, plotly_scale_location_plot, plotly_res_vs_leverage)
# In un chunk diverso per minimizzare cpu-time
# Normalizziamo le covariate
df$f1_z <- scale(df$f1)
df$f2_z <- scale(df$f2)
df$f3_z <- scale(df$f3)
df$f4_z <- scale(df$f4)
df$f5_z <- scale(df$f5)
linModNormalized <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z, data = df)
summary(linModNormalized)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.624 -3.796 0.044 3.243 17.608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.0000 0.1719 238.505 < 2e-16 ***
## f1_z 2.6573 0.3332 7.976 4.99e-15 ***
## f2_z -2.8564 0.6394 -4.467 9.02e-06 ***
## f3_z -17.7693 1.0322 -17.216 < 2e-16 ***
## f4_z 9.1442 1.1191 8.171 1.13e-15 ***
## f5_z -12.6019 0.2738 -46.024 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.985 on 835 degrees of freedom
## Multiple R-squared: 0.8401, Adjusted R-squared: 0.8392
## F-statistic: 877.6 on 5 and 835 DF, p-value: < 2.2e-16
# Installa e carica le librerie necessarie
if (!requireNamespace("ggplot2", quietly = TRUE)) {
install.packages("ggplot2")
}
if (!requireNamespace("gridExtra", quietly = TRUE)) {
install.packages("gridExtra")
}
library(ggplot2)
library(gridExtra)
# Funzione per creare il plot di Residuals vs Fitted (Residui vs. Valori Previsti)
res_vs_fitted <- ggplot(linModNormalized, aes(.fitted, .resid)) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Residuals vs Fitted", x = "Fitted Values", y = "Residuals") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot QQ dei Residui
qq_plot <- ggplot(linModNormalized, aes(sample = .stdresid)) +
stat_qq(color = "skyblue", size = 2) +
stat_qq_line(color = "red", size = 1) +
labs(title = "QQ Plot dei Residui", x = "Quantiles teorici", y = "Quantiles osservati") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot Scale-Location
scale_location_plot <- ggplot(linModNormalized, aes(.fitted, sqrt(abs(.stdresid)))) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Scale-Location Plot", x = "Fitted Values", y = "sqrt(|Standardized Residuals|)") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot Residuals vs Leverage (Residui vs. Leverage)
res_vs_leverage <- ggplot(linModNormalized, aes(.hat, .stdresid)) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Residuals vs Leverage", x = "Leverage", y = "Standardized Residuals") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Visualizza i quattro plot
grid.arrange(res_vs_fitted, qq_plot, scale_location_plot, res_vs_leverage, ncol = 2)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#PLOT INTERATTIVI DA TESTARE
# Installa e carica le librerie necessarie
if (!requireNamespace("ggplot2", quietly = TRUE)) {
install.packages("ggplot2")
}
if (!requireNamespace("gridExtra", quietly = TRUE)) {
install.packages("gridExtra")
}
if (!requireNamespace("plotly", quietly = TRUE)) {
install.packages("plotly")
}
library(ggplot2)
library(gridExtra)
library(plotly)
# Funzione per creare il plot di Residuals vs Fitted (Residui vs. Valori Previsti)
res_vs_fitted <- ggplot(linModNormalized, aes(.fitted, .resid)) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Residuals vs Fitted", x = "Fitted Values", y = "Residuals") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot QQ dei Residui
qq_plot <- ggplot(linModNormalized, aes(sample = .stdresid)) +
stat_qq(color = "skyblue", size = 2) +
stat_qq_line(color = "red", size = 1) +
labs(title = "QQ Plot dei Residui", x = "Quantiles teorici", y = "Quantiles osservati") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot Scale-Location
scale_location_plot <- ggplot(linModNormalized, aes(.fitted, sqrt(abs(.stdresid)))) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Scale-Location Plot", x = "Fitted Values", y = "sqrt(|Standardized Residuals|)") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Funzione per creare il plot Residuals vs Leverage (Residui vs. Leverage)
res_vs_leverage <- ggplot(linModNormalized, aes(.hat, .stdresid)) +
geom_point(color = "skyblue", alpha = 0.7, size = 3) +
geom_smooth(se = FALSE, color = "red", method = "loess", size = 1) +
labs(title = "Residuals vs Leverage", x = "Leverage", y = "Standardized Residuals") +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
# Converti i singoli grafici in oggetti plotly
plotly_res_vs_fitted <- ggplotly(res_vs_fitted)
## `geom_smooth()` using formula = 'y ~ x'
plotly_qq_plot <- ggplotly(qq_plot)
plotly_scale_location_plot <- ggplotly(scale_location_plot)
## `geom_smooth()` using formula = 'y ~ x'
plotly_res_vs_leverage <- ggplotly(res_vs_leverage)
## `geom_smooth()` using formula = 'y ~ x'
# Visualizza i quattro grafici interattivi
subplot(plotly_res_vs_fitted, plotly_qq_plot, plotly_scale_location_plot, plotly_res_vs_leverage)
# TEST SUL MODELLO DI REGRESSIONE LINEARE
#1 Summary
summary (linModNormalized)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.624 -3.796 0.044 3.243 17.608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.0000 0.1719 238.505 < 2e-16 ***
## f1_z 2.6573 0.3332 7.976 4.99e-15 ***
## f2_z -2.8564 0.6394 -4.467 9.02e-06 ***
## f3_z -17.7693 1.0322 -17.216 < 2e-16 ***
## f4_z 9.1442 1.1191 8.171 1.13e-15 ***
## f5_z -12.6019 0.2738 -46.024 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.985 on 835 degrees of freedom
## Multiple R-squared: 0.8401, Adjusted R-squared: 0.8392
## F-statistic: 877.6 on 5 and 835 DF, p-value: < 2.2e-16
#2 R-quadrato e R-quadrato Adattato
summary_linModNormalized <- summary(linModNormalized)
r_squared <- summary_linModNormalized$r.squared
cat("R-squared:", r_squared, "\n")
## R-squared: 0.8401246
n <- length(df$o_oreb)
k <- length(linModNormalized$coefficients) - 1
adjusted_r_squared <- 1 - ((1 - r_squared) * (n - 1) / (n - k - 1))
cat("Adjusted R-squared:", adjusted_r_squared, "\n")
## Adjusted R-squared: 0.8391672
#2 test Shapiro per valutare la normalita' dei residui
shapiro.test(residuals(linModNormalized))
##
## Shapiro-Wilk normality test
##
## data: residuals(linModNormalized)
## W = 0.99508, p-value = 0.008216
#3 test di omoschedasticita' (Breusch-Pagan test) --> risultato suggerisce omoschedasiticita'
bptest(linModNormalized)
##
## studentized Breusch-Pagan test
##
## data: linModNormalized
## BP = 9.6069, df = 5, p-value = 0.08717
#4 test di multicollinearita'
car::vif(linModNormalized)
## f1_z f2_z f3_z f4_z f5_z
## 3.751989 13.820282 36.007839 42.333607 2.534096
Possiamo ora verificare la presenza di outlayers, ovvero dati anomali o fuori norma che possono influenzare in modo molto significativo il modello rispetto al resto dei dati presi in cosiderazione. Provvediamo quindi a rimuoverli e a ricreare il modello.
influencePlot(linModNormalized,5)
## StudRes Hat CookD
## 541 3.1086834 0.002953697 0.004722471
## 769 -0.3942749 0.044140265 0.001197642
## 822 3.6069525 0.027320444 0.060040657
## 957 2.9775866 0.024929163 0.037426310
#rimuovo gli outlayers e ricreo il modello
df_1 <- df[-which(row.names(df) %in% c(822, 957)),]
df_1$f1 <- (df_1$o_oreb)/(df_1$o_fga-df_1$o_fgm)
df_1$f2 <- (df_1$d_dreb)/(df_1$d_fga-df_1$d_fgm)
df_1$f3 <- (df_1$o_oreb + 1.5 * df_1$d_dreb)/(df_1$o_dreb + 2 * df_1$d_oreb)
df_1$f4 <- (df_1$o_oreb - df_1$o_dreb) + 1.5 * (df_1$d_dreb - df_1$d_oreb)
df_1$f5 <- (df_1$d_oreb / df_1$d_to) / (df_1$o_dreb / df_1$o_to)
df_1$f1_z <- scale(df_1$f1)
df_1$f2_z <- scale(df_1$f2)
df_1$f3_z <- scale(df_1$f3)
df_1$f4_z <- scale(df_1$f4)
df_1$f5_z <- scale(df_1$f5)
linModNormalized_1 <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z, data = df_1)
summary(linModNormalized)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.624 -3.796 0.044 3.243 17.608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.0000 0.1719 238.505 < 2e-16 ***
## f1_z 2.6573 0.3332 7.976 4.99e-15 ***
## f2_z -2.8564 0.6394 -4.467 9.02e-06 ***
## f3_z -17.7693 1.0322 -17.216 < 2e-16 ***
## f4_z 9.1442 1.1191 8.171 1.13e-15 ***
## f5_z -12.6019 0.2738 -46.024 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.985 on 835 degrees of freedom
## Multiple R-squared: 0.8401, Adjusted R-squared: 0.8392
## F-statistic: 877.6 on 5 and 835 DF, p-value: < 2.2e-16
summary(linModNormalized_1)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z, data = df_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.8618 -3.7270 -0.0235 3.2212 15.5573
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.0417 0.1700 241.386 < 2e-16 ***
## f1_z 2.4792 0.3322 7.462 2.15e-13 ***
## f2_z -3.1936 0.6314 -5.058 5.21e-07 ***
## f3_z -18.0997 1.0167 -17.803 < 2e-16 ***
## f4_z 9.6926 1.1055 8.767 < 2e-16 ***
## f5_z -12.6741 0.2714 -46.698 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.925 on 833 degrees of freedom
## Multiple R-squared: 0.8432, Adjusted R-squared: 0.8423
## F-statistic: 896.1 on 5 and 833 DF, p-value: < 2.2e-16
# Installa e carica le librerie necessarie
if (!requireNamespace("ggplot2", quietly = TRUE)) {
install.packages("ggplot2")
}
if (!requireNamespace("ggfortify", quietly = TRUE)) {
install.packages("ggfortify")
}
library(ggplot2)
library(ggfortify)
## Warning: il pacchetto 'ggfortify' è stato creato con R versione 4.3.2
# Crea il grafico di influenza
influence_plot <- influencePlot(linModNormalized_1, id = 5, main = "Influence Plot")
## Warning in applyDefaults(id, defaults = list(method = "noteworthy", n = 2, :
## unnamed id arguments, will be ignored
# Personalizza il grafico
#influence_plot +
# theme_minimal() + # Usa uno stile minimale
# theme(plot.title = element_text(size = 16, face = "bold"), # Titolo più grande e in grassetto
# axis.title = element_text(size = 14), # Etichette degli assi più grandi
# axis.text = element_text(size = 12), # Numeri degli assi più grandi
# legend.position="none") # Nascondi la legenda
# #INTERATTIVO DA TESTARE
# # Installa e carica le librerie necessarie
# if (!requireNamespace("car", quietly = TRUE)) {
# install.packages("car")
# }
# if (!requireNamespace("plotly", quietly = TRUE)) {
# install.packages("plotly")
# }
# library(car)
# library(plotly)
#
# # Calcola l'influenza
# influence_values <- influence(linModNormalized_1)
#
# # Estrai i valori influenti per l'osservazione desiderata (es. id=5)
# influence_values_single <- as.matrix(influence_values$infmat[5, ])
#
# # Crea un grafico interattivo
# plot_ly(x = colnames(influence_values_single), y = influence_values_single,
# type = "bar", name = "Influence", marker = list(color = "blue")) %>%
# layout(title = "Influence Plot",
# xaxis = list(title = "Variables"),
# yaxis = list(title = "Influence"))
Nel complesso entrambi i modelli sembrano essere buoni, spiegano entrambi circa l’84% della varianza nella variabile dipendente “won”, anche se la rimozione degli otlayers ha migliorato leggermente i risultati.
Dato che la differenza tra i due modelli e’ molto ridotta (<0.5%), tale da essere trascurabile, abbiamo deciso di utilizzare il modello comprendente gli outlayers.
# Divisione in Test e Train per evitare che il modello fitti troppo bene sui nostri dati
sample <- sample(c(TRUE, FALSE), nrow(df), replace=TRUE, prob=c(0.7, 0.3))
train <- df_reb[sample, ]
test <- df_reb[!sample, ]
# METODO DI REGRESSIONE LASSO
#define response variable
y <- train$won
#define matrix of predictor variables (uso solo poche variabili ma potete farlo con tutte da togliere però la risposta area)
x <- data.matrix(train)
#perform k-fold cross-validation to find optimal lambda value, la cross-validation è un ottimo modo per non overfittare e trovare il miglior modello
cv_model <- cv.glmnet(x, y, alpha = 1)
#find optimal lambda value that minimizes test MSE
best_lambda <- cv_model$lambda.min
best_lambda
## [1] 0.364785
#produce plot of test MSE by lambda value
plot(cv_model)
# Fittiamo il modello con il best lambda (penalizzazione)
best_model <- glmnet(x, y, alpha = 1, lambda = best_lambda)
coef(best_model)
## 8 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 1.1833088
## o_oreb .
## o_dreb .
## o_reb .
## d_oreb .
## d_dreb .
## d_reb .
## won 0.9708495
# SANDBOX TESTING
# QUI SI TESTA TUTTO CIÒ CHE POI VERRÀ AGGGIUNTO NEL DOCUMENTO SOPRA
values <- aggregate(cbind(o_oreb, o_dreb, d_oreb, d_dreb, o_reb, d_reb, won) ~ tmID, data = df, FUN = sum)
# temp <- hist (temp, col = 'steelblue', main = 'caccaculo', xlab = 'balls')
### SANDBOX
teams <- c()
df$reb <- c(df$o_reb + df$d_reb)
summary(df$reb)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5963 6805 7002 7040 7245 8359
for (team in df$tmID)
{
teams <- unique(c(teams, team))
}
teams <- sort(teams)
values <- aggregate(cbind(reb, o_oreb, o_dreb, d_oreb, d_dreb, o_reb, d_reb, won) ~ tmID, data = df, FUN = sum)
x <- seq(-4, 4, by = 0.01)
y <- dnorm(x, mean = 0, sd = 1)
plot (x,y, type="l")
normal <- data.frame(x, y)
df$reb <- df$o_reb + df$d_reb
# ISTOGRAMMA INTERATTIVO
plot_ly(data = df, x = ~reb, type = "histogram",
marker = list(color = 'lightblue', line = list(color = 'black', width = 1)),
nbinsx = 20) %>%
layout(title = 'Istogramma dei Rimbalzi',
xaxis = list(title = 'Rimbalzi'),
yaxis = list(title = 'Frequenza'),
shapes = list(type = "line",
x0 = mean(df$reb), x1 = mean(df$reb),
y0 = 0, y1 = 215,
line = list(color = "red", width = 2, dash = "dash")),
yaxis = list(autorange = TRUE))
#PLOT DI DENSITA CON SOVRAPPOSIZIONE DI UNA NORMALE IDEALE
density_data <- density(df$reb)
mu <- mean(df$reb)
sigma <- sd(df$reb)
normal_data <- dnorm(density_data$x, mean = mu, sd = sigma)
p <- plot_ly(x = density_data$x, y = density_data$y, type = 'scatter', mode = 'lines',
line = list(color = 'blue', width = 2),
name = "Densità rimbalzi totale") %>%
layout(title = "Densità dei Rimbalzi",
xaxis = list(title = "Rimbalzi"),
yaxis = list(title = "Density", autotick = TRUE, autorange = TRUE))
p <- add_trace(p, x = density_data$x, y = normal_data, mode = 'lines',
line = list(color = 'green', width = 2),
fill = "tozeroy", fillcolor = "rgba(0, 255, 0, 0.2)",
name = "Dist normale ideale")
p <- add_trace(p, x = c(mu, mu), y = c(0, max(density_data$y)),
mode = 'lines', line = list(color = 'red', width = 2, dash = 'dash'),
name = "Media")
(p)
# CORRPLOT
data_subset <- df[, c("reb", "o_reb", "d_reb")]
cor_matrix <- cor(data_subset)
rownames(cor_matrix) <- colnames(data_subset)
colnames(cor_matrix) <- colnames(data_subset)
cor_data <- reshape2::melt(cor_matrix)
names(cor_data) <- c("Var1", "Var2", "Corr")
dimnames(cor_matrix) <- list(rownames(cor_matrix), colnames(cor_matrix))
data_for_plotly <- as.data.frame(as.table(cor_matrix))
p <- plot_ly(data = cor_data,
x = ~Var1,
y = ~Var2,
z = ~Corr,
type = "heatmap",
colors = colorRampPalette(c("#4575b4", "#91bfdb", "#e0f3f8", "#fee08b", "#d73027"))(100),
hoverinfo = "x+y+z") %>%
layout(title = 'Correlation Matrix',
xaxis = list(title = "", tickangle = 45, side = "bottom", automargin = TRUE),
yaxis = list(title = "", automargin = TRUE),
autosize = TRUE)
cor_values <- round(as.matrix(cor_matrix), 2) # Round for readability
for (i in seq_len(nrow(cor_matrix))) {
for (j in seq_len(ncol(cor_matrix))) {
p <- p %>% add_annotations(
x = rownames(cor_matrix)[i],
y = colnames(cor_matrix)[j],
text = as.character(cor_values[i, j]),
showarrow = FALSE,
font = list(color = ifelse(cor_values[i, j] < 0.5, "white", "black"))
)
}
}
(p)
#BOXPLOT (1° Versione)
aggregated_data <- df %>%
select(tmID, year, o_reb, d_reb) %>%
group_by(tmID, year) %>%
summarize(
o_reb = mean(o_reb, na.rm = TRUE),
d_reb = mean(d_reb, na.rm = TRUE)
)
## `summarise()` has grouped output by 'tmID'. You can override using the
## `.groups` argument.
reshaped_data <- aggregated_data %>%
pivot_longer(cols = c(o_reb, d_reb), names_to = "stat_type", values_to = "stat") %>%
unite("new_col", tmID, stat_type, sep = "_") %>%
pivot_wider(names_from = new_col, values_from = "stat")
fig <- plot_ly()
for (i in 2:ncol(reshaped_data)) {
current_column_data <- reshaped_data[[i]]
fig <- fig %>% add_trace(y = current_column_data, name = colnames(reshaped_data)[i], type = "box")
}
#BOXPLOT (2° Versione)
# Creare categorie per i rimbalzi
df$reb_categories <- cut(df$reb, breaks = c(-Inf, 5, 10, 15, Inf), labels = c("0-5", "6-10", "11-15", "16+"))
# Convertire df$won in numeri
df$won <- as.numeric(as.character(df$won))
# Creare un boxplot con ggplot2
ggplot(df, aes(x = reb_categories, y = won, fill = reb_categories)) +
geom_boxplot(color = "black", notch = FALSE) +
labs(title = "Distribuzione delle Vittorie in base ai Rimbalzi",
x = "Rimbalzi",
y = "Vittorie") +
theme_minimal()
#TEST ANDERSON-DARLING
ad.test(df$reb)
##
## Anderson-Darling normality test
##
## data: df$reb
## A = 3.6997, p-value = 3.1e-09
#Con un livello di significatività (α) di 0.01 e un p-value molto piccolo (3.1e-09) ottenuto dal test di normalità di Anderson-Darling per i dati della variabile df\(reb, puoi concludere che hai sufficiente evidenza statistica per respingere lipotesi nulla che i dati seguono una distribuzione normale.Con il tuo livello di significatività del 0.01 e il p-value molto piccolo (3.1e-09), il p-value è inferiore al livello di significatività, quindi respingeresti lipotesi nulla. Questo suggerisce che i dati nella variabile df\)reb non seguono una distribuzione normale al livello di significatività del 0.01. In termini più pratici, hai abbastanza evidenza statistica per concludere che la variabile df$reb non segue una distribuzione normale basandoti sui risultati del test di Anderson-Darling.
#TEST KOLMOGOROV SMIRNOV
ks.test(df$reb, "pnorm")
## Warning in ks.test.default(df$reb, "pnorm"): i legami non dovrebbero essere
## presenti per il test di Kolmogorov-Smirnov
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: df$reb
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided
#TEST SHAPIRO WILK
sf.test(df$reb)
##
## Shapiro-Francia normality test
##
## data: df$reb
## W = 0.98016, p-value = 1.758e-08
#In sintesi, il risultato del test di Shapiro-Francia indica che i tuoi dati nella variabile df$reb non seguono una distribuzione normale. Questo è supportato dal valore basso del p-value, il quale suggerisce che la differenza tra la distribuzione dei tuoi dati e una distribuzione normale è statisticamente significativa.
barplot(df$reb, col = c("#1b98e0", "#353436"))
heatmap_df <- subset(values, select = -c(o_oreb, o_dreb, d_oreb, d_dreb, won, reb))
rownames(heatmap_df) <- heatmap_df$tmID
heatmap_df <- heatmap_df[,-1]
heatmaply(
heatmap_df,
colors = viridis(n = 256, option = "magma"),
k_col = 2,
k_row = 4,
)
plot_ly(values, x = ~tmID, y = ~o_reb, type = 'bar', name = 'Rimbalzi offensivi', marker = list(color = '#FFAFA1')) %>%
add_trace(y = ~d_reb, name = 'Rimbalzi difensivi', marker = list(color = '#b2fff8')) %>%
layout(yaxis = list(title = 'Valori'), barmode = 'stack')
#VARIABILI DUMMY
# VARIABILI DUMMY CON TEST ANOVA su ConferenceID
# Controlla i valori univoci nella colonna 'conference'
unique_conferences <- unique(df$confID)
print(unique_conferences)
## [1] "EC" "WC"
# Supponiamo che 'conference' sia la variabile categorica di interesse
df$conference <- as.factor(df$confID)
# Creazione di variabili dummy
dummy_vars <- model.matrix(~ df$confID - 1, data=df) # -1 per evitare la colinearità
# Aggiunta delle variabili dummy al dataframe
df <- cbind(df, dummy_vars)
# Esecuzione di un test ANOVA
# Supponiamo che 'won' sia la variabile dipendente
anova(lm(df$won ~ df$confID, data=df))
## Analysis of Variance Table
##
## Response: df$won
## Df Sum Sq Mean Sq F value Pr(>F)
## df$confID 1 252 251.61 1.6295 0.2021
## Residuals 839 129548 154.41
# si tiene l'ipotesi nulla perchè la variabile correlazione non ha un effetto sulle vittorie.
# Non Posso inserire a modello la variabile
# VARIABILI DUMMY CON TEST ANOVA su DivisionID
# Controlla i valori univoci nella colonna 'conference'
unique_conferences <- unique(df$divID)
print(unique_conferences)
## [1] "CD" "AT" "MW" "PC" "SE" "SW" "NW"
# Supponiamo che 'conference' sia la variabile categorica di interesse
df$conference <- as.factor(df$divID)
# Creazione di variabili dummy
dummy_vars <- model.matrix(~ df$divID - 1, data=df) # -1 per evitare la colinearità
# Aggiunta delle variabili dummy al dataframe
df <- cbind(df, dummy_vars)
# Esecuzione di un test ANOVA
# Supponiamo che 'won' sia la variabile dipendente
anova(lm(df$won ~ df$divID, data=df))
## Analysis of Variance Table
##
## Response: df$won
## Df Sum Sq Mean Sq F value Pr(>F)
## df$divID 6 2723 453.78 2.9782 0.006964 **
## Residuals 834 127077 152.37
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# si rifiuta l'ipotesi perchè la variabile divisione ha un effetto sulle vittorie. Posso inserire a modello la variabile
linModNormalized1 <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z + divID, data = df)
summary(linModNormalized1)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z + divID,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.5087 -3.7310 -0.2418 3.2627 17.7887
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.0693 0.3623 113.356 < 2e-16 ***
## f1_z 2.7683 0.3338 8.294 4.39e-16 ***
## f2_z -2.5867 0.6410 -4.035 5.95e-05 ***
## f3_z -17.4386 1.0717 -16.273 < 2e-16 ***
## f4_z 8.6520 1.1528 7.505 1.58e-13 ***
## f5_z -12.5564 0.2730 -45.998 < 2e-16 ***
## divIDCD -0.3150 0.4995 -0.631 0.5285
## divIDMW -0.6553 0.5318 -1.232 0.2182
## divIDNW -1.7753 1.0766 -1.649 0.0995 .
## divIDPC 0.7862 0.5098 1.542 0.1234
## divIDSE -1.8621 1.0808 -1.723 0.0853 .
## divIDSW 2.1069 1.0993 1.917 0.0556 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.945 on 829 degrees of freedom
## Multiple R-squared: 0.8438, Adjusted R-squared: 0.8418
## F-statistic: 407.2 on 11 and 829 DF, p-value: < 2.2e-16
#EFFETTI INTERAZIONE
linModNormalized1 <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z + divID + f1_z:f2_z + f4_z:divID, data = df)
summary(linModNormalized1)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z + divID +
## f1_z:f2_z + f4_z:divID, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.7588 -3.5284 -0.1167 3.2668 16.8359
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.7756 0.3810 107.009 < 2e-16 ***
## f1_z 2.6552 0.3340 7.950 6.16e-15 ***
## f2_z -2.6926 0.6397 -4.209 2.84e-05 ***
## f3_z -17.6878 1.0989 -16.096 < 2e-16 ***
## f4_z 8.2411 1.2243 6.731 3.15e-11 ***
## f5_z -12.5192 0.2729 -45.872 < 2e-16 ***
## divIDCD -0.2156 0.4998 -0.431 0.66630
## divIDMW -0.4261 0.5348 -0.797 0.42582
## divIDNW -1.7830 1.0962 -1.626 0.10423
## divIDPC 0.8067 0.5097 1.583 0.11384
## divIDSE -1.7857 1.2263 -1.456 0.14575
## divIDSW 2.2319 1.1597 1.925 0.05462 .
## f1_z:f2_z -0.4993 0.1764 -2.831 0.00475 **
## f4_z:divIDCD 1.1625 0.5271 2.205 0.02770 *
## f4_z:divIDMW 1.4944 0.4977 3.003 0.00276 **
## f4_z:divIDNW 1.7033 1.5238 1.118 0.26395
## f4_z:divIDPC 0.4753 0.4939 0.962 0.33614
## f4_z:divIDSE 0.3874 1.2052 0.321 0.74795
## f4_z:divIDSW -0.2418 1.6248 -0.149 0.88176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.911 on 822 degrees of freedom
## Multiple R-squared: 0.8473, Adjusted R-squared: 0.844
## F-statistic: 253.4 on 18 and 822 DF, p-value: < 2.2e-16
# f4 è indicatore che cresce all'aumentare dei rimbalzi ottenuti e diminuisce all'aumentare dei rimbalzi subiti.
# Interazione tra f4 e la divisione non è significativa, contrariamente al fatto che invece interazione tra f1 e f2 è significativo
# AGGIORNAMENTO Del modello lineare togliendo la seconda interazione
linModNormalized2 <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z + divID + f1_z:f2_z, data = df)
#test F, significatività singoli coefficienti, R^2
summary(linModNormalized2)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z + divID +
## f1_z:f2_z, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.3485 -3.6885 -0.2249 3.0529 18.2413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.7526 0.3810 106.955 < 2e-16 ***
## f1_z 2.7181 0.3332 8.158 1.26e-15 ***
## f2_z -2.6976 0.6402 -4.214 2.79e-05 ***
## f3_z -17.8375 1.0789 -16.533 < 2e-16 ***
## f4_z 9.1576 1.1651 7.860 1.20e-14 ***
## f5_z -12.5511 0.2720 -46.136 < 2e-16 ***
## divIDCD -0.1705 0.5009 -0.340 0.73363
## divIDMW -0.4737 0.5345 -0.886 0.37582
## divIDNW -1.5260 1.0772 -1.417 0.15696
## divIDPC 0.8806 0.5094 1.729 0.08422 .
## divIDSE -1.8421 1.0771 -1.710 0.08760 .
## divIDSW 2.1476 1.0956 1.960 0.05032 .
## f1_z:f2_z -0.4492 0.1727 -2.601 0.00946 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.928 on 828 degrees of freedom
## Multiple R-squared: 0.8451, Adjusted R-squared: 0.8429
## F-statistic: 376.4 on 12 and 828 DF, p-value: < 2.2e-16
#Nel complesso, il modello sembra ben adattato ai dati, spiegando una parte significativa della variazione nelle vittorie delle squadre di pallacanestro. La significatività delle variabili derivate (f1_z, f2_z, f3_z, f4_z, f5_z) e delle variabili dummy divID suggerisce che questi fattori sono importanti nella predizione del numero di vittorie. Il termine di interazione f1_z:f2_z è significativo, indicando che c'è un effetto sinergico tra le variabili f1_z e f2_z nella loro influenza sul numero di vittorie.
linModNormalized2_pois = glm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z + divID + f1_z:f2_z, family=poisson(link=log), data = df)
summary(linModNormalized2_pois)
##
## Call:
## glm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z + divID +
## f1_z:f2_z, family = poisson(link = log), data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.673408 0.012235 300.247 < 2e-16 ***
## f1_z 0.056118 0.010780 5.206 1.93e-07 ***
## f2_z -0.088471 0.020809 -4.252 2.12e-05 ***
## f3_z -0.507549 0.036027 -14.088 < 2e-16 ***
## f4_z 0.298790 0.038518 7.757 8.68e-15 ***
## f5_z -0.326690 0.009074 -36.002 < 2e-16 ***
## divIDCD -0.007013 0.015957 -0.440 0.6603
## divIDMW -0.028657 0.017160 -1.670 0.0949 .
## divIDNW -0.036591 0.035179 -1.040 0.2983
## divIDPC 0.008163 0.016123 0.506 0.6127
## divIDSE -0.038462 0.035741 -1.076 0.2819
## divIDSW 0.042715 0.032817 1.302 0.1931
## f1_z:f2_z -0.010771 0.005525 -1.950 0.0512 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3362.04 on 840 degrees of freedom
## Residual deviance: 595.35 on 828 degrees of freedom
## AIC: 5248.2
##
## Number of Fisher Scoring iterations: 4
normal = coefficients(linModNormalized2)
poisson = exp(coefficients(linModNormalized2_pois))
cbind(normal, poisson)
## normal poisson
## (Intercept) 40.7526017 39.3858990
## f1_z 2.7180646 1.0577226
## f2_z -2.6975972 0.9153300
## f3_z -17.8374742 0.6019693
## f4_z 9.1575970 1.3482271
## f5_z -12.5510548 0.7213074
## divIDCD -0.1705107 0.9930112
## divIDMW -0.4736569 0.9717494
## divIDNW -1.5259888 0.9640706
## divIDPC 0.8805990 1.0081959
## divIDSE -1.8421079 0.9622686
## divIDSW 2.1475512 1.0436400
## f1_z:f2_z -0.4492010 0.9892864
#Considerazioni Generali:
#È importante notare che gli effetti delle variabili possono essere interpretati in modo diverso a seconda della distribuzione scelta per il modello. La scelta tra distribuzione normale e di Poisson dipende dalla natura della tua variabile dipendente e dai tuoi obiettivi di modellazione. Nel contesto di modelli di regressione, è sempre buona pratica verificare l'adeguatezza del modello esaminando i residui, eseguendo test diagnostici e valutando la bontà di adattamento. L'interpretazione dei coefficienti dovrebbe essere fatta considerando la scala appropriata per la distribuzione utilizzata (lineare per la normale, logaritmica per la Poisson). Se stai cercando di prevedere il numero di vittorie, la distribuzione di Poisson potrebbe essere più appropriata per variabili conteggio come questa. Tuttavia, è sempre necessario verificare l'adeguatezza del modello ai dati specifici.